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I.  INTRODUCTION/PHYSICAL  DESCRIPTION  OF  THE  PROBLEM 


With  the  rapidly  increasing  development  of  higher  power  lasers, 
and  the  growing  number  of  applicat?ons  of  these  devices,  a number 
of  new  and  transient  heat  transfer  problems  require  solutions. 

One  such  set  of  problems  involves  tnc  absorption  of  energy  from 
a laser  beam  by  the  gas  through  which  the  beam  is  passing.  This 
report  describes  a numerical  technique  by  which  the  various  energy 
transfer  and  gas  motion  phenomena  in  such  an  irradiated  gas 
can  be  modeled.  The  problem  consists  of  a gas  mixture  enclosed 
in  a vessel  of  square  cross-section.  The  beam  passes  through  the 
gas  mixture  along  the  center  of  the  square  section  as  chown  in 
Figure  1.  The  container  is  considered  to  be  well  insulated  and 
sealed.  This  implies  that  the  laser  energy  is  incident  through  a 
sealed  window  or  an  open  hole  in  the  container  which  is  very  small 
compared  to  the  volume  of  the  container.  Thus  the  pressure  within 
the  enclosure  is  unrelieved  and  increases  as  the  gas  is  heated. 

The  beam  diameter  is  small  compared  to  the  cross-section  dimensions 
of  the  enclosure  and  the  amount  of  energy  absorbed  is  relatively 
small . 

The  medium  or  fluid  being  heated  consists  of  a homogeneous  mixture 
of  gases.  The  gas  composition  and  chemical  properties  do  not  change. 

The  energy  absorbed  is  assumed  to  manifest  itself  only  as  a change 
in  the  temperature  of  the  gas  with  no  chemical  change.  For  longer 
wavelength  radiation,  the  absence  of  chemical  change  due  to  photochemical 
effects  is  probably  a reasonable  assumption.  However,  as  the  gas 
mixture  becomes  hotter,  some  chemical  changes  will  begin  to  occur. 

The  temperatures  at  vhich  these  chemical  effects  become  significant 


depends  upon  the  reactivity  of  the  gases  in  the  mixture.  Therefore, 
this  analysis  is  not  applicable  where  temperatures  are  high  enough 
to  cause  chemical  changes  in  the  gas  mixture.  The  mixture  con- 
sidered in  this  analysis  is  a normal  hydrocarbon  vapor  in  air. 
Absorption  coefficients  for  these  vapors  are  extracted  from 
Reference  1,  and  the  air  is  assumed  to  be  completely  transparent 
to  the  laser  radiation.  (The  ratio  of  the  absorption  coefficients 
of  air  to  pentane  is  10  ^.)  These  absorption  coefficient  apply  to 
10.6  micron  radiation  as  would  be  produced  by  a carbon  dioxide 
laser,  and  are  assumed  to  be  directly  proportional  to  the  gas 
density  (i.e.,  number  of  molecules  in  beam  path).  As  a first  order 
approximation,  the  absorption  coefficient  is  assumed  to  be  independent 
of  temperature.  Other  gas  mixtures  can  be  included  in  this  analysis 
if  sufficient  thermal  properties  and  absorption  coefficients  at  the 
specified  wavelength  are  known  for  each  constituent.  Similarly,  other 
wavelength  laser  radiation  can  be  considered. 

This  analysis  of  the  problem  consists  of  the  construction  of  i 
computer  program  which  models  the  heat  transfer  effects  occurring  in 
the  gas  (i.e.,  convection,  conduction,  and  absorption  of  the  laser 
radiation).  Re-radiation  from  the  gas  is  not  considered.  This  should 
not  be  a problem  at  the  relatively  low  temperatures  to  which  th  gases 
are  heated.  Since  only  a very  3mall  amount  of  energy  is  removed 
from  the  beam  as  it  passes  through  the  gas,  the  gas  mixture  temperature 
gradients  and  incident  energy  gradients  along  tho  beam  axis  are  very 

small.  Therefore,  the  problem  is  assumed  to  be  two-dimensional. 
Furthermore,  symmetry  about  a vertical  plane  through  che  center  of 


the  b^am  is  assumed. 


3 


Two  different  types  of  laser  beams,  characterized  by  different 
spatial  energy  distributions,  are  considered.  The  first  of  these 
is  the  "flat-top"  or  "square"  beam  which  has  a constant  power  density 
versus  radial  distance  from  the  beam  center  distribution,  as  shown  In 
Figure  2.  The  second  type  of  beam  considered  has  a Gaussian  power 
density  distribution  truncated  at  the  defined  radius.  This  is  also 
shown  in  Figure  2.  The  laser  output  power  is  assumed  to  be  constant 
with  time,  i.e.,  continuous  wave  (CW)  output. 

II.  MATHEMATICAL  DESCRIPTION  OF  THE  PROBLEM 

The  coordinate  system  utilized  in  this  analysis  is  shown  in 
Figure  3.  Terms  used  u this  report  are  defined  in  Appendix  1. 

A.  ASSUMPTIONS: 

The  problem  is  assumed  to  be  two-dimensional  and  symmetric 
about  a vertical  plane  through  the  laser  beam  axis.  The  beam  passes 
through  the  center  of  the  square  enclosure.  Only  the  right  side  of 
the  cross-section  of  the  enclosure  will  be  considered  as  shown  in 
Figure  3.  The  enclosure  is  assumed  to  be  insulated  and  the  system 
is  therefore  adiabatic,  with  the  exception  of  the  volumetric  heating 
of  the  gas  by  absorption  of  energy  from  the  laser  beam.  The  density 
of  the  gas  is  assumed  to  remain  constant  except  in  the  absorption 
coefficient  calculation  and  in  the  bouyancy  term  in  the  vertical 
equation  of  motion.  This  assumption  is  discussed  in  Appendix  3. 

The  heating  occurs  in  a sealed  enclosure.  Thus,  there  is 
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no  pressure  relief.  The  average  gas  density  (p)  in  the  enclosure 
remains  constant  and  equal  to  the  initial  gas  density  (Po)  . The 
specific  heat,  thermal  conductivity,  viscosity,  and  absorption 
coefficient  of  the  gas  mixture  remain  constant.  The  gas  mixture  is 
assumed  to  act  as  a perfect  gas,  and  the  viscous  energy  dissipation 
terms  in  the  energy  equation  are  neglected. 

E.  INITIAL  AND  BOUNDARY  CONDITIONS 


The  gas  mixture  is  initially  at  rest  and  at  initial  tem- 
perature (0O)  , pressure  (P0) , and  density  (po).  At  time  t = 0 
the  laser  beam  is  turned  on.  The  boundary  conditions  are: 

1.  Zero  velocity  at  walls 

at  x = 0 and  x=L:  u = v = 0 

at  y = L/2 : u = v * 0 

2.  Zero  temperature  gradient  at  walls  (adiabatic) 

at  x = 0 and  x = L:  — = 0 

9x 


at  y = L/j 


3.  Zero  mass  flow  and  zero  gradients  of  temperature 
vertical  velocity  across  line  of  symmetry 


and 


at  y = 0:  V = 0,  fy-  = 0,  fy  = 0 


C.  EQUATIONS 


1.  x component  of  motion  equation: 
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2.  y component  of  motion  equation: 
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3.  Continuity  equation: 
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4.  Energy  equation: 

30  . 30  . 30  k ,320  . 30,  , q 

at  * “ + v F J1^ 


9x 


ay 


III.  METHOD  OF  SOLUTION: 

A.  DERIVATION  OF  EQUATIONS 

The  method  of  deriving  the  required  equations  is  similar 
to  that  of  Reference  2.  Differentiating  the  ; component  of  the 


equation  of  motion  with  respect  to  y,  yields: 


2 2 2 
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(Eq.  1) 


Similarly,  differentiating  the  y component  of  the  equation  of 

motion  with  respect  to  x,  yields: 

2 2 2 2 
~ .9  - . 9u  9v  , 9 v , 9 v , 9v  . 

P St Ox  + 9x  9x  + U 2 + 9x9y  + V 9x9y) 

3X 

2 2 2 
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3x3y  3x  3x2  3y2 

Subtracting  equation  2 from  equation  1:  ^ 
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(Eq.  2) 


Since  = 0 (continuity  equation)  , we  see  that  the  second  and 

3 x 3y 

fourth  terms  on  the  left  side  of  the  equation  cancel,  yielding: 
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(Eq.  3) 
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This  equation  of  motion  and  the  energy  equation  will  now  be  non- 
dimensional  via  the  following  relationships: 


Coordinates : X = — 
r 

Velocities:  U = — 

v 

0 

Temperature:  T = — - 1 


Time  : 


t = 


tv 


Y 

r 

V 


Density:  <4  = — ™ 4- 
Po  P 


Heating  rate : Q = — — 


2 

Pc  vF 


Applying  these  relationships  to  equation  3 yields: 


p [—  — (—  M 
1 2 3t  < 2 3Y 
r r 


r 3 2 
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v 32V 
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v 3 U 


+ UL\  _ I 1_  /V_  32V  , V 32v, 

3 „ 2J  r 3X  1 2 + T — 
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Simplifying  and  using  - = v 
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(Eq.  4) 


Applying  the  non-dimensionalizing  relationships  to  the  energy 
equation  produces: 


|I  + V v io_  3T  V Q_0  3T  = k_  A 0O  32T. 

r2  r M h r23x2+7^ 

+ h «“9„ 

r 

Simplifying  again  using  the  definition  of  kinematic  viscosity,  and 

using  the  Prandtl  number  (Np^  = Ui.  ) > yields 

II  + n II  . w II  J^T  32T.  , 

3t  U 3X  + V 3Y  N ^ 2 + 2 + ^ 

PR  3X  3Y 
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(Eq.  5) 


i . .1.— . 


Finally,  applying  the  non-dimensionalizing  relationships  to  the 

initial  and  boundary  conditions  produces  the  following: 

at  t = 0 : U = V = 0,  T = 0 

T 

at  X = 0 and  X = — : U = V = 0,  ^-=0 

r dX 


L ST 

31  Y = 2r  : U - V = 0 , ^ = 0 


_3T 

3Y 


at  Y = 0:  = 0,-  - = 0,  V = 0 


The  stream  function  (tjj)  and  vorticity  (?)  are  now  introduced: 


U-f  (Eq.  6) 


- 1 (E<>-  7> 


S = ~ A = . 3_|  . 3^  (Eq.  8) 

3X  3Y 

Applying  these  stream  function  and  vorticity  equations  to  the 
motion  equation  (equation  4)  yields: 


+ U + V = Er  A4.  4.  /p  g\ 

3t  U 3X  V 3Y  2 3Y  V C (Eq’  9) 

v 

Equations  5 through  9 describe  the  problem,  and  the  boundary  and 
initial  conditions  have  now  become  the  following: 

at  i=0:  U = V = ? = 0,  T - 0 

T 

at  X = 0 and  X = — : U=V=0,  -^=0 

r 3X 

at  Y = ^ : U = V = 0,  ||  = 0 

3Y  ’ 3Y  3X  4 


B.  SOLUTION  OF  THE  EQUATIONS 
Laser  Heating  Function  Q 

Equations  5-9  can  now  be  solved  via  a series  of  finite 

10 


difference  approximations  advancing  in  time.  The  driving  function  is 

the  non-dimensional  volumetric  heating  rate  term  (Q)  in  the  energy 

equation.  For  small  absorption  coefficients,  the  actual  volumetric 

heating  rate  (q)  is  approximately  the  product  of  the  incident  power  den- 

2 

sity,  I(X,Y)  and  the  absorption  coefficient  a(X,Y,x).  I(X,Y)  (watts/cm  ) 

varies  with  location  in  accordance  with  the  formulas  in  Figure  2. 

The  absorption  coefficient  (o , cm  is  assumed  to  be  dependent  only 

on  the  local  gas  density,  P(X,Y,x),  which  varies  with  location  and 

time.  Thus,  the  actual  volumetric  heating  rate  must  be  calculated 

at  each  grid  point  (X,  Y)  in  the  beam  radius  and  at  each  time  t . The 

local  density  distribution  at  time  t is  used  to  calculate  the  heating 

rate  distribution  at  time  (t  + At  ) . 

q (X,Y , T + At)  - I(X,Y)  <x(X,Y,x) 

q is  then  non-dimensionalized  via  the  previously  mentioned 

formula:  2 

Q (X, Y,  t + At)  = 3-^ — q (X,Y,  t + At)  (Eq.  10) 

Pcv6 

o 

Energy  Equation 

The  next  step  requires  the  solution  of  the  energy  equation 
to  find  the  new  distribution  of  local  temperatures  at  time  (t  + At  ) . 
Using  standard  finite  difference  approximations  with  a forward 
differenced  time  derivative,  centered-dif f erence  spatial  derivatives, 
and  a Taylor  series  approximation  ror  the  second  order  terms,  the 
energy  equation  (equation  5)  can  be  expressed  as: 

^ [T(X,Y,  t + Ax)  - T(X,Y,t)  ] + U (X,Y,t)  [T(X  + AX.Y.t) 

- T(X  - AX.Y.t)]  + V (X,Y,t)  ~ [T(X,Y  + AY,t)  - T(X,Y  - AY,t)] 

- [T(X  + AX,Y,t)  + T(X  - AX.Y.t)  + T(X,Y  + AY,t) 

NPR  (AX)2 


11 


+ T(X,Y  - AY,t)  - 4T(X,Y, t) ] + Q(X,Y,  x + At) 


where  AX  = AY 


Since  the  distributions  of  all  the  variables  are  known  at.  time  x,  and 
Q(X,Y,  x+  At)  can  be  determined  from  equation  10,  we  may  solve  for 
T(X,Y,t  + At  ) . 

T(X,Y,  t + At)  = — ~ [T(X  + AX.Y.t)  + T(X  - AX,Y,i) 

PR  (AX) 

+ T(X,Y  + AY,t)  + T(X, Y - AY , t ) - 4T(X,Y,x)]  + AtQ(X,Y,t  + At) 


At 


- U(X.Y.t)  ~ [T(X  + AX.Y.t)  - T(X  - AX.Y.t)] 


At 


-V  (X.Y.x)  — [T(X,Y  + AY.t)  - T(X,Y  - AY.x)]  + T(X,Y,t)  (Eq.  11) 


Since  the  first  derivative  of  T at  the  walls  and  across  the  line 
of  symmetry  is  zero,  the  values  of  T along  these  boundaries  may  be 
defined  by  using  a standard  "mirror  image"  technique,  [e.g., 
at  T(X,o , t ) : T(X,Y  + A Y,  t)  = T(X,Y  - A Y, t ) ] . 

Density  Function 

An  average  pressure  term,  P(x  ) , is  determined  at  timer  + Ax 
in  the  following  manner: 

Assuming  the  gas  mixture  to  be  a perfect  gas. 


P(x  + Ax)  - P(t)  = ~ [0 (t  + At)  - 0(x)] 


(Eq.  12) 


MV 


where  0 is  the  mass  average  temperature  of  the  total  mass  of  the  gas 
(m)  in  the  enclosure  having  a constant  volume  V.  This  implies  that 
the  pressure  is  homogeneous  within  the  enclosure  or  that  pressure 
equalization  occurs  much  more  rapidly  than  temperature  equalization. 


Furthermore,  the  total  heat  input  during  any  particular  time 
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— 


■ 


interval,  At,  is  At  //  q(X,Y, T+AT>dXdY,  which  can  also  be  expressed 
Y X 

At  Z I __  ..  . , \ 


q(X,Y, t + Ax) 


This  total  heat  input  can  be  related  to  the  mass  average  temperature 


rise  in  the  following  manner: 


Y x q(x»Y»  T + At)  = [0(t  + At)  - 0 (t) ] 


(Eq.  13) 


Substituting  for  m [0  (t  + At)  - 9(t)]  in  equations  12  and  13,  and 


simplifying,  yields: 


P(t  At)  - P(t)  = — l l q(X ,Y , t + At)  (Eq.  14) 

MVc  Y X 


which  can  be  solved  for  P(t  + At) . 


. . . v P (X, Y,  t + At)  P (t  + At)  9o 

*(X-Y.  t + At)  - - p 0(X>Y>  T + At) 


a »0 

from  T = — - 1,  it  is  clear  that  — = [T(X,Y,  t + At)  + 1] 
9o  9 


substituting  into  the  equation  for  tj>,  there  follows: 


d>(X,Y,  t + At)  = 


P(t  + At) 


P [T(X,Y , T + At)  + 1] 
o 


(Eq.  15) 


Vorticity  Equation 


Equations  10,  11,  14,  and  15  can  be  evaluated  in  succession 


to  yield  values  of  Q,  T,  P,  and  $ , respectively,  at  the  advanced  time 


x + At-  The  vorticity  function  (X,Y,t)  is  now  advanced  to  time 


t + At  via  the  motion  equation.  The  motion  equation  (eq.  9)  can  be 


expressed  numerically  in  the  following  manner: 

7-  U(X,Y,  t + At)  - C(X,Y,t)]  + U(X,Y,t)  777  U(X  + AX.Y.t) 


- C(X  - AX.Y.t)]  + V(X,Y , t)  [C(X,Y  + AY,t)  - C(X,Y  - AY,t)] 


= —7  [<KX,Y  + AY, t)  - $ (X,Y  - AY , t ) ] + I C (X  + AX.Y.t) 


+ C(X  - AX.Y.t)  + C (X,Y  + AY, t)  + t(X,Y  - AY,t)  - 4c(X,Y,t)] 
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—ft-.  - - - 


This  can  be  solved  for  5 (X,Y,  t+  At). 


?(X,Y,  t + At)  = ?(X,Y,t)  - U(X,Y,t)[c(X  + AX.Y.t) 

- C(X  - AX,Y,t)]  - |^-  V(X,Y, x)  [ £(X,  Y + AY,t)  - s(X,Y  - AY,t)] 


+ S~  2AY  [4>(X,Y  + AY’  T + At)  “ *(X’Y  " AY*  T + At>  + 
v (AX) 

U(X+AX,Y,t)  + c(X-AX,Y,t)  + ?(X,Y+AY,t)  + s(X,Y-AY,t)  - 4c(X,Y,t)] 

(Eq.  16) 

Equation  16  can  be  solved  to  find  the  values  of  vorticity  at  all 
interior  points.  Along  the  line  of  symmetry,  the  vorticity  function 
is  zero.  A parabolic  extrapolation  of  the  vorticity  values  approach- 
ing any  of  the  three  walls  is  used  to  establish  vorticity  values 
along  these  walls.  This  extrapolation  technique  is  described  below 
and  was  used  in  Reference  3. 

Consider  the  5 values  approaching  the  vertical  wall  for  a constant 


X value,  as  shown  below: 


Y Y Y 
vt-3  w-2  w-1 


where  the  separation  between  each  of  the  Y values  is  AY.  A parabola 
can  be  constructed  to  intersect  the  three  points  (5  , Y ) , (r  , ,v  _i  ) » 

\v  W X W JL 

(^w_2*  yw-2^  • This  parabola  can  be  defined  by  the  following  equation: 

(Y  - Y ) (Y  - Y ) £ (Y  ) 

C,  (Y) 5-* ^ 

2 (ay  r 

(Y  - Y ,)  (Y  - Y ) C(Y  .)  (Y  - Y ) (Y  - Y ) ?(Y  ) 
w— 2 w w— 1 w— 2 w— 1 w 

(AY)2  2 (AY) 2 


-133 


Differentiation  of  this  equation  yields: 

(2Y  - Y - Y ) ?(Y  ,) 

w-1  w w-2 


C(Y) 


2(AY)‘ 


(2Y  - Y - Y ) C(Y  .)  (2Y  - Y _ - Y .)  £(Y  ) 

w-2  w w-1  w-2  w-1  w 

(AY)2  2 (AY) 2 

A 

This  equation  an  now  be  used  to  determine  ? (Y)  at  Y = \_2' 


‘ (V2> 


(2Y  , - Y , - Y ) ;(Y  ,) 

w-2  w-1  w w-2 

.2 


2(AY)‘ 


(2Y  „ - Y , - Y ) ?(Y  .)  (2Y  - Y - Y ) ?(Y  ) 

- w-2  w-2  w s w-1  , w-2  w-2  w-1 w_ 

2 + 9 

2(AY)“ 


(AY)‘ 


Utilizing  a centered  difference  approximation  for  5 - i.e.,  £(Y  ^) 

“ 2(|y)"  ^^w-1^  “ 5(Yw_3)land  remembering  that 


Y - 
w 

Yw-2 

2AY 

Y - 
w 

Vi  = 

AY 

Y - Y 
w-1  w-2 


AY 


the  equation  is  simplified  to  produce: 

2AY  I^(Yw-l)  " ^(Yw-3)]  " 2AY  C(YW-2)  + AY  ^(Yw-l)  ~ 2AY 


Solving  this  equation  for  s( Y ),  one  obtains: 


5<V  - ScCY^j)  - 3C(Yw_2)  + C(Yv_3)  (Eq.  17) 

Thus,  the  value  of  the  vorticity  is  determined  along  the  vertical 
wall.  A similar  procedure  is  used  for  the  lower  and  upper  walls. 

Stream  Function 

Now  that  the  values  of  vorticity  are  known  at  time  t + At, 

the  stream  function  (i|>)  can  be  determined.  This  is  done  by  numerically 

solving  the  defining  equation  for  vorticity  (eq.  8),  for  if).  Using 
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truncated  Taylor  series  approximations  for  the  second  order  terms, 
equation  8 can  be  expressed  as: 


C(X,Y,t)  = - 


[||/(X  + AX.Y.t)  + 4>(X  - AX,  Y,t)  - 2i|y(X,Y,T)] 


[*(X,Y  + AY, t)  + i|»(X,Y  - AY,t ) - 2<KX,Y,r)] 


Although  the  above  equation  can  be  solved  via  the  explicit 


"relaxation"  technique,  the  number  of  calculations  required  before 


suitable  convergence  occurs  is  excessive.  Thus,  an  alternating 


direction  implicit  (ADI)  method  was  used  to  solve  the  equation  for 


, This  ADI  method  is  used  to  find  a solution  to  the  equation 

li  = ijk.  + ijfe.  + , 

-A.,  O T T *» 


2 2 
3X  3Y 


where  the  term  y is  analogous  to  a time  term.  The  objective  is  to 


perform  iterations  until  the  term  ^ approaches  zero  (analogous  to  a 


steady-state  solution),  thereby  providing  a solution  to  the  vorticity 


2 

equation  r,  — V 1/1 


The  ADI  computation  proceeds  in  steps  of  -^Ay  . The  partial 


differential  equation  is  written  in  finite  difference  form,  looking 


forward  with  respect  to  y and  implicit  in  the  X direction  during  the 


first  step  of  -Ay  . A second  equation,  implicit  in  the  Y direction, 


is  then  used  for  the  second  st^p  of  y Ay  . The  procedure  is  repeated 


until  the  change  in  ip  across  a pair  of  steps  of  ^Ay  is  sufficiently 


close  to  zero,  thereby  yielding  a steady-state  solution  with  respect 
to  Y . We  now  introduce  # and  ip^Y)  > representing  the  value 


of  $ at  timer  , the  value  of  after  solution  of  the  first  ADI  equation. 


and  the  value  of  ^ after  solution  of  the  second  ADI  equation,  respectively. 


Thus  the  two  ADI  equations  are  as  shown  below. 


Implicit  in  X direction: 


- 


($Ay)  (o)  GAy)  GAy)  GAy) 

»(X,Y)  - »(X,Y)  = »(X+AX,Y)  - 2»(X,Y)  + »(X-&X,Y) 

*AY  (AX)2 

(o)  (o)  (o) 

+ iL^Y^JjL.(X?Y).  f ».(X,Y-AY)  * ,(X>Y>  t+At> 


(Eq.  18) 


Implicit  in  Y direction: 

(Ay)  GAy)  GAy)  GAy)  GAy) 

»(X,Y)  - »(X,Y)  _ »(X+AX.Y)  - 2»(X,Y)  + iHX-AX.Y) 

2 


(Ay)  (Ay)  (Ay) 

+ »(X,Y+  Y)  - 2fr  (X.Y)  + »(X, Y-AY) 
2 


- + C(X,Y,  t + At) 
Gay)  ^ 


(Eq.  19) 


Equation  18  can  be  solved  for  i>  , and  the  results  used  in 

/ .v \ (Ay)  (o) 

equation  19  to  determine  . If  *(X,Y)  - iKX.Y)  is  sufficiently 

(Ay)  (Ay)  (o) 

small  then  <KX,Y,  t + At)  is  set  equal  to  iJi(X,Y).  If  ip(X,Y)  “ 't'(X.Y) 

(o)  (Ay) 

is  too  great,  then  'I'(X.Y)  is  set  equal  to  iKX.Y)  and  the  calculation 
is  repeated.  An  arbitrary  value  of  Ay  may  be  chosen;  however, an 
optimum  value  exists  at  which  computation  time  will  be  minimized.  The 
method  by  which  equations  18  and  19  are  solved  will  now  be  described. 

Multiplying  equation  18  by  ($0^  (remembering  that  AX  = AY)  and 
collecting  like  terms,  we  obtain; 

GAy)  r . Y \ 2 GAY)  GAy) 


GAy)  ,ax,2  GAy)  (iAy) 

i|»(X+AX9Y)  - 2 (1  + ) <KX,7)  + *(X-AX,Y) 


(o)  (o) 


- (AX)2  q(X,Y,  t + At)  -<KX,Y+AY)  - 2G|~i-  - 1)  *(X,Y)-*(X,Y-AY) 

(Eq.  20) 

This  equation  will  be  solved  over  a column  of  grid  points  (Y  = constant) 

Thus,  the  series  of  equations  is  of  the  following  form: 

(iAy)  (iAy)  (iAy) 

ip(xi+AX,Y)  - B«j.(Xi,Y)  + <)(Xi-AX,Y)  = c(X1)  i = 1,M  (Eq.  21) 


where  B * 2 (1  + ~^r)  , a known ^constant 


C * right  hand  side  of  eq.  20,  a known  constant  for  each 

and  M * ^ + ] , the  number  of  grid  points  in  the  column. 

Thus,  there  is  a set  of  M simultaneous  linear  equations  with  M 
(iAy) 

unknown  values  of  >HX^,Y)  . Since  the  velocities  at  the  walls  are  zero, 
’J'  * constant  along  the  walls.  Furthermore,  since  the  line  Y * 0 is 
a line  of  symmetry,”  iKX.Y,  t)*i|>(X,-Y,t)  , which  implies  tha*  t|»(X,0,t)  * 
0.  Since  the  line  of  symmetry  intersects  the  walls  which  have  ij/  * con- 
stant, it  follows  that  the  value  of  along  the  walls  is  also  zero. 

(*Ay) 

Thus,  the  values  of  ^(X^Y)  in  equation  21  are  already  known  for  i * 1 
and  i * M.  The  set  of  equations  represented  by  equation  21  may  be 
written: 

(iAY)  (*Ay) 

0 ~ B<|»  (X2,Y)  + *(X3,Y)  - c2 

(iAy)  (iAy)  ($Ay) 

<KX2Y)  - Bi>  (X3,Y)  + <kx4,y)  - c3 


(iAy)  (jAy)  (iAy) 

+ ♦(XH_1,Y).  cm_2 

(iAy)  (i Ay ) 

*(XM-2’Y'  - B* <Vl’Y)  + 0 ‘ Vl 

which  can  be  expressed  in  matrix  form  as: 


“ ™ 

(iAy) 

- *“ 

• B 1 0 0 ...  0 

*(x2,y) 

c2 

1 -B  1 0 ...  0 

(i  A y ) 

*(x3,y) 

C3 

0 1 -B  1 ...  0 

0 

m 

C4 

• 

0 

• 

0 

| * • • • • • 

* • • 

0 

0 

o 1 -B  1 

(i  Ay ) 

0 

0 0 1 -B 

<,<xm-i,y) 

CM-] 

h 
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The  tri-diagonal  square  matrix  can  then  be  inverted  and  values  of 
(iAy) 

(X  ,Y)  can  be  obtained  for  a given  column.  Equation  19  may  be  solved 

('»!  (iAy) 

in  exactly  the  same  manner  to  obtain  ij/(X,Y)  in  terms  of  (X * Y j ) 

along  rows  (X  * constant,  j M 1 to  ~~  + 1).  The  procedure  is 

repeated  until  a steady-state  solution  occurs,  at  which  point  the  new 

(Ay) 

stream  function  sji  (X,Y,  t + At  ) is  set  equal  to  <KX,Y). 

Velocity  Function 

The  final  step  in  the  solution  of  this  set  of  equations  for 
a single  time  step  advancement  At  , is  the  determination  of  the  non- 
dimensional  velocities  from  the  stream  function  \J>  . This  is  accomplished 
via  solution  of  equations  6 and  7.  A centered-difference  approximation 
is  used  and  the  equat  .ons  to  be  solved  are  as  shown  below: 

U(X,Y,  t + At)  = ~ U<X»Y+AY.  T+At)  - *(X,Y-AY,  t+At)]  (Eq  - 22) 

V(X,Y,  T + At)  = ~ [iKX+AX,Y,  t+At)  - iKX-AX,Y,  t+At)]  (Eq.  23) 

At  this  point,  all  parameters  have  been  advanced  to  the  new  time 
t + At,  and  the  cycle  is  ready  to  begin  again.  A flow  chart  of  this 
procedure  is  shown  in  Figure  U. 

Conservative  and  Upwind  Differencing 
Several  runs  were  made  with  51  x 26  and  101  x 51  grids 
representing  four  inch  and  e.;ht  inch  cross-sectional  areas  respectively 
Qx  =Ay  * .08").  Beam  diameters  were  varied  from  one  to  three  inches, 
and  time  steps  ( At)  of  from  .0001  to  .001  seconds  were  used.  Although 
the  overall  solutions  appeared  reasonable,  apparent  discretization 
errors  did  occur,  producing  unreasonably  cool  regions  in  the  tempera- 
ture arrays.  The  most  obvious  remedy  to  this  problem  would  be  the 

lntroductirn  of  still  finer  grids  and/or  smaller  time  steps.  However, 
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because  the  memory  storage  required  for  all  these  arrays  and  the  com-. 

(' 

puter  time  required  was  already  quite  large,  another  approach  was 
taken.  This  approach  consisted  of  the  addition  of  conservative  and 
upwind  differencing  numerical  approximations  to  the  energy  and  motion 
equations . 

The  conservative  difference  approximation  is  derived  from  the 
product  rule  for  differentiation,  and  incorporates  the  condition  of 
mass  conservation  in  the  differentiation  process.  To  demonstrate  its 


application,  let  us  consider  the  second  and  third  terms  ir  the  non- 


9T  9T 

dimensionalized  energy  equation  (eq.  5) , U gx  + ^ gY  • Utilizing  the 


product  rule  for  differentiation. 


U + V = — — (UT)  - T — + — (VT)  - T — 
3X  3Y  3X  K * 3X  3Y  V ' 9Y 


3U  av 

and  since  3v  + Tv  0 from  continuity, 


.,31  . v 3T  _ 9 (UT)  3 (VT) 

3X  3Y  3X  3Y 


which  may  be  expressed  numerically  as: 


2^  f U (X+AX , Y , t ) T(X+AX,Y,t)  - U(X-AX,Y,t)  T(X-AX,Y,x)] 


+ 2 AY  fv(X’Y+AY»T)  T(X,Y+AY,t)  - V(X,Y-AY,r)  T(X,Y-AY,t)] 


The  above  expression  can  now  be  substituted  into  the  energy 


equation  (eq.  11).  The  second  and  third  terms  on  the  right  hand  side 


H a c, 

of  the  equation  of  motion  (eq.  16),  representing  U ^ + V jy 


are  similarly  converted  to  the  conservative  numerical  form. 

The  upwind  differencing  technique  is  also  applied  to  the  same 
terms  of  the  energy  and  motion  equations.  This  technique  replaces 
the  conservative,  centered-dif f eronce  , numerical  approximation  with  a 
backward  looking  (with  rnspect  to  velocity),  conservative,  numerical 


approximation.  That  is,  the  derivative  is  taken  in  the  direction 





"from  which  the  wind  is  blowing".  For  example,  the  term 
is  approximated  by 


8(UT) 

ax 


^ [U(X,Y,t)  T(X,Y,t)  - U(X-AX,Y,t)  T(X-AX,Y,t) ] 
if  U (x,y,t)  is  positive,  and 

[U(X+AX,Y,t)  T(X+AX,Y,t)  - U(X,Y,t)  T(X,Y,t)] 
if  U(x,y,r)  is  negative. 

The  application  of  the  upwind,  conservative , numerical  approxima 
tions  to  the  appropriate  terms  of  the  energy  and  motion  equations, 
resulted  in  smoother  and  physically  reasonable  temperature  and  vorticity 
arrays.  The  unreasonably  cool  zones  in  the  temperature  distribution 
were  eliminated  and  the  lines  of  constant  vorticity  became  smoother. 

C.  GAS  MIXTURE  PROPERTIES 

One  portion  of  the  input  d;  '.a  requirements  is  the  properties 
of  the  gas  mixture.  The  properties  of  each  component  (up  to  ten  com- 
ponents) are  input  to  the  program  and  are  used  to  calculate  properties 
of  the  mixture.  The  following  are  the  properties  of  each  component 
which  are  input : 

Density  (p)  at  initial  temperature  (T0)  and  pressure  (P0) 
Absorption  coefficient  (a)  per  psi  per  volume  percent  at 
initial  density  (pG) 

Thermal  conductivity  (k) 

Molecular  weight  (M) 

Molar  concentration  (mole  fraction)  (n) 

Kinematic  Viscosity  (v) 

Specific  Heat  (c) 

The  method  by  which  these  properties  are  used  to  calculate 

the  gas  mixture  properties  will  now  be  described.  In  the  following 

22 


discussion,  the  subscripts  i and  j refer  to  the  individual  gas  mixture 
components'  properties,  while  m refers  to  the  mixture  properties. 


where  P is  viscosity 


N 


m 


R0 


n M = p 
i i o 


i = 1 


N = no.  of  components 


N 


M 

m 


N 


i=  1 


niMi 


a = 100  P 'S  n.a. 

o , i i 


N 


i = 1 


1 

' 5~  2 ci"iMi 

m i = 1 

The  viscosity  of  the  mixture  is  determined  via  the  semi- 
empirical  formula  of  Wilke  as  described  in  Reference  4.  That  is; 


N 

2 
i = 1 


Vi 


where 

A 


ij 


n1Aij 

j --  1 ] 3 


and  of  course, 
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The  thermal  conductivity  of  the  mixture  is  determined  by  a similar 
method,  attributed  to  Mason  and  Saxena,  and  also  described  in  Reference 
4.  That  is; 
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where  A^  is  the  same  as  in  the  viscosity  equation. 


IV.  APPLICATION  EXAMPLES 

This  section  consists  of  a sample  set  of  input  Cata  and  the  results 
of  four  runs  made  with  slight  variations  in  this  data.  Appendix  ?.  is 


a listing  of  the  main  program  and  the  four  primary  subroutines. 

Instructions  for  the  use  of  this  program  are  given  in  Appendix  4 
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A.  SAMPLE  PROBLEM  INFUT/OUTPUT: 

Input  data  can  be  divided  into  two  subsets:  problem  parameters 

and  gas  mixture  constituents'  properties.  These  are  listed  in  Figure 
5,  along  with  the  values  assigned  for  the  first  sample  problem.  Since 
the  thermal  conductivity,  viscosity,  and  specific  heat  remain  fixed 
throughout  the  program,  despite  local  temperature  increases  of  a 
several  hundreds  of  degrees  Rankine,  the  values  of  those  three  proper- 
ties shown  in  Figure  5 are  actual  and/or  estimated  values  at  approxi- 
mately 672°R  (100°C)  . The  universal  gas  constant  and  the  acceleration 
of  gravity  are  constant  values  contained  within  the  program.  Only  very 
minor  changes  are  required  to  vary  the  acceleration  of  gravity  or  the 
spatial  power  distribution  as  is  done  in  two  of  the  examples. 

Output  data  is  presented  both  graphically  and  in  tabular 
form.  The  tabular  output  data  presented  in  the  following  example 
problems  are  the  properties  of  the  gas  mixture,  the  actual  beam  power 
and  area  used  in  the  calculations  (which  are  different  from  the  input 
conditions  due  to  approximations  to  the  heated  area  and  due  to  trunca- 
tion of  Gaussian  power  distribution  beams) , the  average  power  density 

9 

or  intensity  (watts/cm  ) , the  pressure  in  the  enclosure,  and  the  rate 
at  which  the  beam  is  heating  a one-inch  thick  section  of  tie  gas 
mixture.  The  graphical  output  data  consists  of  contour  plots  of 
temperature,  stream  function  (streamlines)  and  vorticity  within  the 
enclosure,  and  of  vertical  velocity  along  three  horizontal  lines  of 
constant  X.  All  of  the  above  information  can  be  obtained  at  any  time 
after  irradiation  has  begmi. 

R . SAMPLE  PROBLEM  NO . 1 : 

This  problem  is  described  by  the  input  data  values  listed  in 


PROBLEM  PARAMETERS: 


Initial  Temperature  (deg.  rankine) 

Laser  Beam  Power  (watts) 

Beam  Radius  (inches) 

Time  Interval  - At  (sec) 

Initial  Pressure  (psia) 

Number  of  Components  in  the  Gas  Mixture 
Grid  Spacing  - Ax  or  Ay 

Enclosure  Size  - length  of  side  of  square 
cross-section  (inches) 

Beam  Duration  (sec) 


520. 

63675. 


1.5 


.001 

14.7 

2 


.08 


8.00 


.8 


GAS  MIXTURE  CONSTITUENTS'  PROPERTIES: 
Property 

Name 

Density  (lbm/cu.  ft.) 

Absorption  Coefficient  (cm  ^ per  psi 
per  volume  percent) 

Thermal  Conductivity  (btu/hr.  ft  °R) 
Molecular  Weight 
Concentration  (mole  fraction) 
Kinematic  Viscosity  (sq.  ft/sec) 
Specific  Heat  (btu/lbm  °R) 


Gas  1 
Pentane 
.183 

-5 

1.39  x 10 
.0096 
72.15 
.03 

3.1  x 10-5 
.409 


Gas  2 

Air 

.0735 

0 

.017 

28.97 

.97 

2.18  x 10 
.24 


FIGURE  5:  INPUT  PARAMETERS  & SAMPLE  PROBLEM  VALUES 
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Figure  5.  In  addition,  the  standard  acceleration  of  gravity  and  a 
Gaussian  power  distribution  within  the  beam  are  used.  The  calculated 
properties  of  the  gas  mixture  are  listed  below: 

Density  = .07974  lbm/cu  ft 

Absorption  Coefficient  = .000613  cm-^ 

Thermal  Conductivity  = .01628  btu/ (hr  ft  °R) 

Molecular  Weight  = 30.265 

Kinematic  Viscosity  = .000189  sq  ft/sec 

Specific  Heat  = .252  btu/(lbm  °R) 

The  program  assumes  heating  to  occur  at  each  grid  point  which 
lies  within  the  defined  beam  radius.  Thus  the  heated  area  only 
approximates  a semi-circle.  This  heated  area  is  shown  in  the  graphical 
output.  Also,  the  Gaussian  beam  is  truncated  at  the  defined  radius. 

Thus,  approximately  13%  of  the  Gaussian  laser  beam  power  is  omitted 
from  the  calculation.  As  a result  of  these  approximations  and 
truncations,  the  actual  power  used  was  57,586.8  watts,  the  actual 
heated  area  was  7.283  sq  inches,  and  the  resulting  power  density 
(intensity),  averaged  over  the  entire  heated  area  was  1225.86  watts/sq  cm. 

Figures  6 through  15  are  the  graphical  output  for  this  problem. 
The  roughly  semi-circular  area  drawn  on  these  plots  represents  the 
heated  area.  The  contours  of  each  of  the  variables  plotted  are 
identified  by  symbols.  The  line  identified  by  the  octagon  represents 
the  value  of  the  function  that  lies  ten  percent  of  the  difference 
between  the  maximum  value  and  the  minimum  value  of  that  function, 
above  the  minimum  value.  For  example,  after  .151  seconds  of 
irradiation,  the  element  with  the  highest  temperature  has  a temperature 
of  785. 2°R.  Since  the  lowest  temperature  of  any  element  at  .151 


seconds  is  520°R  (the  initial  temperature) , the  temperature  that  lies 
ten  percent  of  the  way  from  the  lowest  value  to  the  highest  is  546.52°R. 
An  isotherm  of  546.52°R  is  therefore  drawn  and  is  identified  by  an 
octagon.  Similarly,  the  thirty  percent,  fifty  percent,  seventy  per- 
cent, and  ninety  percent  isotherms  are  identiiied  by  a triangle,  a 
plus  sign,  an  x,  and  a diamond,  respectively.  Streamlines  and  lines 
of  constant  vorticity  are  similarly  labelled.  The  plot  of  upward 
Velocity  is  actually  a three-dimensional  drawing  relating  the  non- 
dimensionalized  upward  velocity  at  all  Y values  along  three  lines  of 
constant  X.  The  three  lines  precisely  quarter  the  half  cross-sectional 
area.  Coordinates  and  the  approximate  beam  boundary  are  shewn  in  the 
first  of  these  plots  to  assist  in  visualization. 

At  .151  seconds  (Figures  6 and  7) , the  isotherms  appear  to  be 
nearly  concentric  about  the  centerpoint  of  the  line  of  symmetry,  as 
would  be  expected  with  a Gaussian  power  beam.  However,  there  is  a 
slight  upward  displacement  of  these  isotherms  due  to  the  convective 
motion  of  the  gas  which  is  now  beginning  to  influence  the  pattern. 

The  plots  of  upward  velocity  and  stream  function  demonstrate  this  con- 
vective motion.  Note  that  the  highest  value  lines  of  constant  stream 
function  and  vorticity  are  approximately  centered  about  a point  near 
the  outermost  horizontal  point  of  the  heated  region. 

At  .294  seconds  (Figures  8 and  9),  the  effects  of  convection 
begin  to  appear.  The  regions  of  highest  temperature,  stream 
function,  and  vorticity  have  been  displaced  upward.  The  upward  veloci- 
ties have  increased  significantly  within  and  above  the  heated  region. 

At  this  time  the  maximum  temperature  of  892°R  was  attained.  This  high 
temperature  existed  for  approximately  .02  seconds  and  then  began  to 
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FIGURE  9:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .294  SECONDS  OF  PROBLEM  1 
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FIGURE  11:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .451  SECONDS  OF  PROBLEM  1 
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FIGURE  14:  STREAM  FUNCTION  AND  TEMPERATURE 

AT  .801  SECONDS  OF  PROBLEM  1 
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decrease.  As  with  all  the  other  example  problems,  the  maximum  temper- 
atures always  occur  near  the  intersection  of  the  upper  portion  of  the 
heated  area  with  the  li  .e  of  symmetry. 

At  .451  seconds  (Figures  10  and  11) , the  convective  forces 
have  significantly  affected  the  isotherms.  The  hottest  region  now  lies 
completely  outside  the  heated  area,  and  the  maximum  temperature  has 
decreased  to  850. 7°R.  Most,  and  possibly  all,  of  the  beam  area  is  now 
at  a lower  temperature  than  it  was  at  .294  seconds.  The  flow  is  con- 
tinuing to  accelerate  and  the  maximum  upward  velocity  is  27  inches  per 
second.  This  maximum  velocity  is  occurring  at  a point  on  the  line  of 
symmetry  above  the  heated  region.  Horizontal  motion  is  beginning  to 
alter  the  flow  significantly.  The  lower  portions  of  the  isotherms 
have  become  indented  and  the  vorticity  plot  has  changed  greatly.  The 
central  portion  of  the  streamline  and  constant  vorticity  plots  have 
continued  their  upward  movement. 

At  .601  seconds  (Figures  12  and  13) , the  isothe-ms  have  assumed 
a characteristic  "mushroom  cloud"  shape  and  there  are  two  regions  that 
have  temperatures  greater  than  the  90%  isotherm.  However,  the  maximum 
temperature  has  aecreased  to  770°R.  The  regions  of  highest  stream 
function  and  vorticity  have  moved  upward  and  away  from  the  heated  area. 
The  flow  is  obviously  being  affected  by  the  upper  wall.  The  maximum 
upward  velocity  is  still  increasing,  but  more  slowly,  and  is  now  30.5 
inches /second . 

At  .801  seconds  (Figures  14  and  15),  the  program  is  terrai-^ 
nated.  The  gas  which  was  initially  heated  has  cooled  significantly 
and  been  pushed  to  the  upper  portion  of  the  vertical  wall.  The  hottest 
region  now  lies  within  and  directly  above  the  central  portion  of  the 


beam  area,  along  the  line  of  symmetry.  The  maximum  upward  velocity 
has  decreased  to  about  27  inches/second  and  the  upward  velocity  pro- 
files have  been  significantly  altered.  Apparently,  the  gas  mixture 
is  now  passing  through  the  heated  region  more  slowly.,  as  the  tempera- 
tures in  the  beam  area  have  now  risen  slightly.  The  maximum 
temperature  is  now  792. 5 °R. 

Figure  16  is  a plot  of  the  pressure  rise  within  the  enclosure, 
and  the  heating  rate  per  inch  of  beam  path  through  the  enclosure,  as 
a function  of  time.  The  variation  in  the  heating  rate  is  due  to 
changes  in  the  average  density  of  the  gas  mixture  in  the  beam  path. 

The  absorption  coefficient  is  assumed  to  be  directly  proportional  to 
density.  As  the  density  in  the  heated  region  decreases  due  to  heating, 
the  absorption  coefficient  and  the  heating  rate  also  decrease  propor- 
tionately. Thus,  the  time  at  which  the  heating  rate  is  lowest 
(approximately  .25  seconds)  is  also  the  time  at  which  the  density  in 
the  heated  region  is  lowest.  Since  the  density  increase  due  to  chang- 
ing pressure  is  relatively  small,  this  is  also  the  time  at  which  the 
mass  average  temperature  in  the  beam  path  is  highest. 

C.  SAMPLE  PROBLEM  NO.  2 

This  problem  is  identical  to  problem  number  1,  except  that 
the  laser  beam  spatial  power  distribution  has  been  changed  to  a con- 
stant ("flattop  beam").  Thus,  a heated  element  on  the  edge  of  the 
beam  is  being  irradiated  at  the  same  rate  as  an  element  near  the  center. 
The  input  power  has  been  reduced  so  that  the  total  incident  power  is 
the  same  as  in  problem  number  1.  The  beam  area  is  also  the  same.  The 
program  utilizes  a total  power  of  57603.1  watts  over  an  area  of  7.283 
square  inches,  yielding  an  average  power  density  of  1226.2  watts  per 


HEATING  RATE  (BTU/5EC) 


PRESSURE  RISE  (PSI) 


square  centimeter.  This  is  less  than  .4  watts  per  sq  cm  more  than  the 
average  power  density  of  the  Gaussian  beam  usea  in  the  previous  sample 
problem.  No  other  input  properties  have  changed. 

Although  the  only  change  between  sample  problems  1 and  2 is 
the  manner  in  which  the  power  is  distributed  throughout  the  beam  area, 
the  results  are  significantly  different.  At  time  .151  seconds 
(Figures  17  and  18) , the  isotherms  are  closely  grouped  near  the  edge 
of  the  heated  area.  Higher  thermal,  upward  velocity  and  vorticity 
gradients  exist  near  the  edge  of  the  heated  area  than  at  the  same  time 
in  problem  number  1.  Also,  the  maximum  temperature  rise  is  only  about 
half  that  which  had  occurred  at  this  time  in  problem  1.  The  maximum 
upward  velocity  is  similarly  lower  and  the  maximum  vorticity  is  nearly 
twice  as  great. 

At  .301  seconds  (Figures  19  and  20) , the  effects  of  convection 
are  more  pronounced.  The  ratios  of  the  maximum  temperature,  upward 
velocity,  and  vorticity  at  this  time  to  the  same  values  at  the  same 
time  in  problem  1 have  begun  to  approach  unity;  however,  the  differences 
are  still  quite  large.  For  example,  the  ratio  of  the  maximum  tempera- 
ture in  problem  2 to  the  maximum  temperature  in  problem  1 has  increased 
from  .51  at  .151  seconds  to  .66  at  .301  seconds. 

At  .393  seconds  (Figures  21  and  22)  the  maximum  temperature 
value  of  822°R  occurs.  Notice  that  problem  2 has  resulted  in  a maxi- 
mum temperature  that  is  70°  lower,  and  occurs  0.1  seconds  later  than 
the  maximum  temperature  of  problem  1.  As  in  sample  problem  number  1, 
this  maximum  temperature  value  has  occurred  at  the  top  of  the  heated 

region  near  the  line  of  symmetry. 

41 


As  the  heating  continues,  the  isotherms  assume  the  mushroom 
shaped  pattern  observed  in  problem  1,  and  the  pattern  is  then  simi- 
larly disturbed  by  the  presence  of  the  upper  and  the  vertical  walls. 

The  hottest  region  lies  outside  the  beam  path  for  a longer  time  than 
in  problem  1.  At  .601  seconds  the  maximum  temperatures  and  vorticities 
ai.e  nearly  the  same  as  in  problem  1,  and  the  maximum  velocity  is 
slightly  smellier.  The  plots  of  isotherms,  streamlines,  constant 
vorticities,  and  upward  velocities  at  .601  seconds  and  .801  seconds 
(Figures  23-26)  look  quite  similar  to  the  plots  of  problem  1 at  those 
times.  The  heating  rate  plot  (Figure  27)  is  slightly  different  from 
that  of  problem  1.  The  heating  rate  values  are  slightly  higher  for 
the  flattop  beam,  and  the  minimum  occurs  at  a later  time.  The  net 
effect  is  that  slightly  more  heat  is  added  with  the  flattop  beam, 
despite  the  fact  that  a lower  maximum  temperature  was  obtained,  and 
the  pressure  rise  is  correspondingly  slightly  higher. 

D.  SAMPLE  PROBLEM  NO.  3 

This  problem  is  identical  to  problem  1 with  the  exception  that 
the  acceleration  of  gravity  (g)  has  been  reduced  to  16.1  ft/sec^. 

Since  the  convective  force  term  in  the  equation  of  motion  is  highly 
dependent  on  the  value  of  g,  it  would  be  expected  that  this  change 
would  reduce  the  effects  of  convection.  The  beam  spatial  power  dis- 
tribution for  this  sample  problem  is  the  same  Gaussian  power  distribu- 
tion as  was  used  in  example  1. 

At  .151  seconds  (Figures  28  and  29)  it  is  interesting  to  com- 
pare the  plots  for  this  example  problem  with  those  of  problem  1 
(Figures  6 and  7) . The  isotherm  values  and  shapes  and  the  lines  of 
constant  vorticity  and  stream  function  are  nearly  identical.  However, 


the  values  of  vorticity,  stream  function,  and  upward  velocity  are 
almost  exactly  half  those  of  problem  1.  This  indicates  that  the 
effects  of  convection  are  insignificant  up  to  this  time. 

By  .301  seconds  (Figures  30  and  31)  the  maximum  temperature 
has  reached  945 °R  and  is  still  rising,  whereas  in  problem  1,  the  maxi- 
mum temperature  attained  a maximum  value  of  892°R  at  .293  seconds  and 
began  decreasing.  Throughout  this  low  gravity  problem  the  upward 
velocities  were  lower  thar.  their  corresponding  values  in  the  standard 
gravity  problem.  However,  by  .8  seconds  the  values  were  much  closer 
than  at  early  times. 

At  .384  seconds  (Figures  32  and  33)  the  maximum  temperature  of 
976c  ’R  was  attained.  Thus  the  effect  of  reducing  gravity  to  one-half 
of  its  standard  value  was  to  produce  a 22%  greater  temperature  rise 
(obviously  due  to  the  lessened  convective  cooling  effect)  at  a time 
.09  seconds  later.  By  .601  seconds  (Figures  34  and  35)  the  plots 
appear  almost  identical  to  the  plots  from  problem  1 at  .451  seconds 
(Figures  10  and  11).  However,  the  values  of  these  contours  are  sub- 
stantially different.  As  expected,  the  motion  of  the  gas  has  been 
slowed,  thereby  allowing  the  hot  gas  to  remain  in  the  beam  path  for  a 
longer  time  and  producing  a higher  local  temperature.  Since  the  gas 
within  the  beam  path  is  at  a higher  temperature  (e.g.,  lower  density) 
than  in  problem  1,  and  since  the  absorption  coefficient  is  linearly 
proportional  to  density,  it  is  to  be  expected  that  the  heating  rate 
would  be  lower  for  problem  3.  And,  if  Figure  38  is  compared  with 
Figure  16,  we  find  that  the  heating  rate  is  indeed  lower  for  pioblem  3. 
Thus,  the  lower  heating  rate  has  produced  a higher  temperature. 
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FIGURE  20:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .301  SECONDS  OF  PROBLQt  2 
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FIGURE  21:  STREAM  FUNCTION  AND  TEMPERATURE 

AT  .394  SECONDS  OF  PROBLEM  2 
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FIGURE  22:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .394  SECONDS  OF  PROBLEM  2 
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FIGURE  23:  STREAM  FUNCTION  AND  TEMPERATURE 

AT  .601  SECONDS  OF  PROBLEM  2 
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FIGURE  24:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .601  SECONDS  OF  PROBLEM  2 


Octagon 

Triangle 

+ 

X 

Diamond 
Pine  Tree 


51 


DC 

>— 

</> 


Stream  Function 
Mon-dimensional) 

95.3 

285.9 

476.5 

667.1 

857.7 


Temperature  (°R) 

541.1 

583.2 

625.3 

667.4 

709.5 


Symbol 

Octagon 

Triangle 

+ 

X 

Diamond 


FIGURE  25:  STREAM  FUNCTION  AND  TEMPERATURE 

AT  .801  SECONDS  OF  PROBLEM  2 
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FIGURE  29:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .151  SECONDS  OF  PROBLEM  3 
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FIGURE  33:  VORTICITY  AND  UPWARD  VELOCITY 
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FIGURE  36:  STREAM  FUNCTION  AND  TEMPERATURE 

AT  .801  SECONDS  OF  PROBLEM  3 
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E.  SAMPLE  PROBLEM  NO.  4 


This  sample  problem  is  the  same  as  sample  problem  2,  except 
that  the  beam  radius  has  been  reduced  to  0.5  inches.  The  beam  power 
has  also  been  reduced  in  order  to  maintain  the  same  intensity  as  in 
problem  2,  and  the  beam  power  distribution  is  "flattop".  The  actual 
power  and  total  beam  area  used  in  these  calculations  were  6782.8  watts 
and  .85°  sq.  inches,  respectively.  This  resulted  in  an  intensity  of 
1226.2  watts/sq.  cm. 

At  .151  seconds  (Figures  39  and  40)  the  maximum  temperature 
(654°R)  is  nearly  the  same  as  the  maximum  temperature  of  problem  2 at 
.151  seconds.  A very  simple  calculation  of  the  temperature  rise  of 
the  gas  being  irradiated  was  made.  This  calculation  assumed  no  heat 
loss,  a heating  time  of  .151  seconds,  a heating  rate  which  is  based 
upon  an  average  of  the  heating  rates  at  0,  .05,  .1,  and  .15  seconds, 
and  the  gas  density  at  587°R  (one-half  the  anticipated  rise) . This 
simple  calculation  yielded  a temperature  rise  within  five  degrees  of 
that  calculated  by  the  program.  Thus,  it  appears  that  very  little 
heat  loss  from  the  hottest  region  of  the  heated  area  has  occurred  at 
.15  seconds. 

At  .244  seconds  (Figures  41  and  42)  the  maximum  temperature 
of  712°R  was  attained.  As  in  all  examples,  this  maximum  value  occurred 
near  the  intersection  of  the  upper  portion  of  the  beam  with  the  line 
of  symmetry.  Reducing  the  size  of  the  beam  while  maintaining  intensity 
resulted  in  a decrease  of  approximately  110°R  in  the  maximum  tempera- 
ture attained.  Also,  this  maximum  temperature  occurred  much  sooner 
with  the  small  beam.  For  approximately  the  first  .15-. 20  seconds,  the 
maximum  temperature  values  in  pr^leras  2 and  4 are  nearly  the  same. 
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FIGURE  44:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .451  SECONDS  OF  PROBLEM  4 
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FIGURE  45:  STREAM  FUNCTION  AND  TEMPERATURE 
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FIGURE  46:  VORTICITY  AND  UPWARD  VELOCITY 

AT  .601  SECONDS  OF  PROBLEM  4 
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FIGURE  47.  STREAM  FUNCTION  AND  TEMPERATURE 
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However,  after  .244  seconds  the  maximum  temperatures  for  the  problem 
with  the  large  beam  continue  increasing  for  an  additional  .15  seconds, 
while  the  maximum  temperatures  Cor  the  small  beam  problem  decrease. 

The  upward  velocities  are  very  low  and  slowly  developed  for  the  small 
beam. 

At  .451  seconds  (Figures  43  and  44)  the  characteristic  mush- 
room cloud  shape  is  forming.  The  temperatures  in  the  beam  area  have 
decreased  by  a considerable  amoun..  Alter  this  time,  the  temperatures 
in  the  beam  area  gradually  increase  at  a steadily  decreasing  rate. 

The  ratios  of  the  total  beam  power  and  area  for  problems  2 and  4 are 
the  same,  .1178.  The  ratio  of  the  initial  heating  rates  is  also  .1178, 
and  the  final  heating  rate  ratio,  which  is  approximately  steady-state, 
is  also  .1178. 

V.  CONCLUSIONS,  APPLICATIONS,  AND  RECOMMENDATIONS  FOR  IMPROVEMENT 
A.  CONCLUSIONS 

The  computer  model  desciibed  in  this  report  produces  results 
which  appear  reasonable  and  reliable  for  a Lwo-di mens ion al  heating 
model.  However,  experimental  verificaticn  has  not  been  accomplished. 
Due  to  inaccuracies  and  estimations  in  some  of  the  input  parameters 
and  assumptions,  exact  agreement  with  experimental  results  cannot  be 
expected.  Some  of  these  inaccuracies  and  estimations  will  be  discussed 
latar.  However,  except  for  ti.e  effects  of  chemical  changes/reactions 
in  the  gas  mixture,  the  model  should  produce  at  least  a preliminary 
estimation. 

Based  upon  the  results  of  the  sample  problems  presented  in 


this  report,  and  other  similar  but  unreporce  <.uns  that  were  made. 


it  appears  that  the  effects  of  convection  are  minimal  during  an 
initial  heating  period.  The  duration  of  this  initial  heating  period 
is  dependent  upon  such  parameters  as  the  beam  area  aid  intensity. 

During  this  heating  period  the  primary  mechanism  by  which  the 
absorbed  energy  is  removed  from  the  heated  region  is  conduction. 

For  intensities  and  beam  sizes  such  as  the  ones  considered  in  the 
sample  problems,  conductive  heat  loss  is  negligible  for  the  gas  in 
the  central  portion  of  the  heated  area.  However,  once  sufficient 
energy  has  been  deposited  and  sufficient  time  has  passed,  the  heated 
gas  begins  to  rise  rapidly  enough  to  begin  to  remove  itself  from  the 
beam  area.  It  appears  that  the  highest  temperatures  occur  just  before 
-he  initially  heated  gas  is  swept  from  the  heated  area.  By  this  time, 
convection  is  the  dominant  mode  of  heat  removal.  Prior  to  this  time 
of  highest  temperatures  the  gas  has  been  becoming  increasingly  hotter, 
and  therefore  less  dense.  As  this  occurs,  the  heating  rate  correspond- 
ingly decreases,  due  to  the  decrease  in  the  number  of  potentially 
absorbing  molecules  in  the  beam  path.  The  heating  rate  is  linearly 
proportional  to  the  average  density  of  the  gas  in  the  beam  diameter. 

As  the  convection  motion  of  the  gas  begins  to  remove  the  hot  gas  from 
the  upper  portions  of  the  beam  and  supplies  cooler  gas  from  below  the 
heated  area,  the  average  density  in  the  beam  area  increases.  This 
results  in  an  increase  in  the  heating  rate.  Eventually,  a near 
steady-state  condition  within  the  beam  area  is  establisheu.  Decreasing 
the  acceleration  of  gravity  retards  this  convective  motion  which  then 
allows  the  initially  heated  gas  to  remain  within  the  beam  area  for  a 
longer  time.  Thir  , of  course,  results  in  higher  teraperotures  in  the 
heated  region.  Incre  sing  the  viscosity  of  the  gas  mixture  should 


have  a similar  effect.  The  Gaussian  power  distribution  results  in 
higher  temperatures  near  the  center  of  the  beam,  as  would  be  expected. 
And  finally,  decreasing  the  beam  size  (at  constant  intensity)  results 
in  less  convective  motion  and  increases  the  importance  of  conduction. 
However,  due  to  the  small  beam  size,  much  less  convective  motion  is 
required  to  remove  the  hot  gas  from  the  beam  path.  The  net  effect  of 
decreasing  beam  size  at  constant  intensity  is  to  reduce  the  temperature. 

B.  APPLICATIONS 

This  model  can  be  used  to  determine  how  best  to  maximize  or 
minimize  the  heating  rate,  total  heat  input,  maximum  temperature  and 
the  time  required  to  attain  it,  pressure  rise,  and  reduction  in  tha 
laser  beam  power  as  it  is  transmitted  through  the  gas.  These  studies 
may  in  turn  be  utilized  to  investigate  such  parameters  as  ignition  of 
gases  by  laser  radiation  and  propagation  of  high  energy  laser  beams 
through  the  atmosphere,  within  lasing  gases,  and/or  in  absorbing  gases. 
The  gas  mixture  evaluated  in  this  problem  (3%  pentane  in  air)  is  a 
highly  flammable  mixture.  At  some  temperature  in  excess  of  9bO°R,  this 
mixture  would  be  ignited.  The  exact  temperature  at  which  auto-ignition 
would  occur  is  unknown.  Perhaps  this  model,  coupled  with  experimental 
results,  could  be  used  to  determine  this  temperature.  However, 
any  such  endeavor  would  require  a careful  examination  of  the 
potential  chemical  changes  which  may  be  occurring  as  the  ignition  tem- 
perature is  approached.  Also,  the  potential  occurrence  of  photochemi- 
cal effects  must  be  investigated.  As  a part  of  any  application  of 
this  model,  the  following  variables  should  be  investigated. 

1.  Initial  pressure  and  temperature. 

2.  Beam  parameters  such  as  size,  spatial  power  distribution, 
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and  power.  With  very  slight  modifications,  the  beam  shape 
or  temporal  power  distribution  (e.g.,  pulsed  laser)  can  also 
be  investigated. 

3.  The  acceleration  of  gravity.  This  could  be  a variable  for 
moving  systems  such  as  aircraft. 

4.  Various  gas  mixture  constituents,  as  long  as  suitable 
physical  property  data  can  be  obtained. 

5.  Size  of  the  enclosure.  This  variable  is  probably  not  very 
important  for  short  irradiation  times,  and  the  results 
should  be  applicable  to  non-enclosed  (constant  pressure) 
situations . 

6.  Initial  motion  of  the  gas.  At  the  start  of  the  irradiation, 
it  may  be  possible  (via  modification  of  the  program)  to  give 
the  gas  an  initial  upward  or  downward  velocity.  The  initial 
horizontal  velocity  must  be  zero  to  be  consistent  with 

the  original  assumption  of  symmetry  across  the  verti- 
cal centerline.  This  capability  may  be  useful  for  flowing 
systems . 

7.  Location  of  the  laser  beam  within  the  enclosure.  The 
analysis  described  in  this  report  assumes  that  the  laser 
beam  axis  is  coincident  with  the  axis  of  the  enclosure, 
and  assumes  that  the  enclosure  has  a square  cross-section. 

If  the  laser  beam  is  displaced  vertically,  the  program  can 
be  easily  modified  to  account  for  this.  If  the  laser  beam 
is  displaced  horizontally,  the  assumption  of  symmetry  about 
a central  vertical  plane  is  destroyed,  and  major  changes 

in  the  program  will  be  required.  It  is  anticipated  that 


moving  the  beam  closer  to  the  walls  will  reduce  convective 
motion  and  thereby  affect  the  heating  and  temperature  dis- 
tribution of  the  gas  mixture.  The  effects  of  the  walls 
should  increase  as  the  beam  is  brought  closer  to  the  walls 
and  should  also  increase  with  irradiation  time. 

C.  RECOMMENDATIONS  FOR  IMPROVEMENT 

Potential  thermally  induced  chemical  effects  and  photochemical 
effects  must  be  considered.  These  are  very  briefly  discussed  in 
Section  I.  If  the  absorbed  energy  does  not  quickly  manifest  itself  as 
an  increase  in  the  translational  motion  (temperatures)  of  all  the  mole- 
cules in  the  beam  area,  then  these  effects  may  be  significant.  Also, 
if  the  absorbed  energy  causes  some  chemical  change  in  the  absorbing 
molecules,  or  increases  their  reactivity  independent  of  temperature, 
then  this  must  also  be  considered.  Finally,  as  the  gas  becomes  heated, 
certain  chemical  reactions  may  begin  to  occur  which  release  energy  and 
change  the  chemical  composition  and  physical  properties  of  the  gas 
mixture.  When  this  is  of  concern,  it  may  be  possible  to  couple  some 
chemical  kinetic  model  with  this  heat  transfer  model. 

Another  improvement  which  would  be  of  value  involves  a more 
accurate  approximation  of  the  absorption  coefficient.  For  some  gases, 
the  absorption  coefficient  may  be  a strong  functicn  of  the  total 
absolute  pressure,  the  partial  pressure  of  each  constituent  in  the 
mixture,  and  the  temperature.  These  variations  of  the  absorption 
coefficient  are  in  addition  to  the  dependence  on  density.  Only  the 
effect  of  density  on  absorption  coefficient  is  evaluated  in  this  model. 
Finally,  variations  of  thermal  conductivity,  viscosity,  and  specific 
heat  with  temperature  could  be  incorporated  into  this  model. 
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APPENDIX  1 


LIST  OF  TERMS 


specific  heat  at  constant  pressure  (Btu/lbm°R) 

2 

acceleration  of  gravity  (ft/sec  ) 

2 

local  beam  intensity  or  power  density  (watts/cm  ) 


thermal  conductivity 


width  of  enclosure 


mass  of  gas  in  the  enclosure 
molecular  weight 


mole  fraction 


Prandtl  number 


number  of  components  in  gas  mixture 
pressure  (psi) 

volumetric  heating  rate  (Btu/sec)  of  element 

2 

non-dimensioual  volumetric  heating  rate  = — — 

pc  o 

beam  radius 

universal  gas  constant 

0 

non-dimensional  temperature  = - 1 

y 

O 

rime  after  irradiation  has  begun 
velocity  in  upward  direction 

non-dimensional  upward  density  = — » 

v 3Y 

horizontal  velocity 

vr 

non-dimensional  horizontal  velocity  = — = ~ 


non-dimensional  vertical  coordinate  = - 

r 

non-dimensional  horizontal  coordinate  = ^ 

r 

absorption  coefficient  (cm 

2 

dimensionless  vorticity  * -V  iJj 
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9 * temperature 

9 - mass  average  temperature  of  gas  in  the  enclosure 

v * kinematic  viscosity 

u = viscosity 

P * density  (lbm/ft3) 

P = average  density  within  the  enclosure 

PD  - initial  density  within  the  enclosure  (pQ=  p) 

t = dimensionless  time  = — 

2 

r 

y = dummy  variable,  analogous  to  time 

$ * dimensionless  density  = — — 

, Po 

V * dimensionless  stream  function 


Subscript  o refers  to  initial  conditions 
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APPENDIX  2 


- PROGRAM  LISTING 
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APPENDIX  3 


GAS  DENSITY  VARIATION  ASSUMPTION 

In  the  equations  of  motion,  the  gas  density  is  assumed  to  be 
constant  except  in  the  bouyancy  force  term.  This  assumption  simpli- 
fies the  solution  of  the  equations  and  produces  only  minimal  error 
if  the  variations  in  gas  velocities  and  temperatures  are  small. 
However,  for  some  of  the  problems  considered  here,  local  absolute 
temperature  values  nearly  doubled.  This,  of  course,  results  in  a 
fifty  percent  decrease  in  density.  This  change  in  density  due  to 
heating  is  considered  in  the  calculation  to  determine  the  amount  of 
energy  absorbed  from  tne  laser  beam,  but  is  not  taken  into  account 
on  the  left  side  of  the  equations  of  motion.  Since  the  magnitude  of 
the  error  introduced  by  this  assumption  will  vary  spatially  with 
temperature,  its  effects  should  be  most  significant  near  the  heated 
region.  Unfortunately,  this  is  also  the  area  of  greatest  interest. 
It  is  recommended  that  future  attempts  to  solve  this  problem  incor- 
porate variable  density  throughout. 


104 


J - 1-. 


APPENDIX  4 


PROGRAM  USERS'  INSTRUCTIONS 


Figure  5 lists  the  input  data,  and  the  appropriate  units  for  this 


data,  which  must  be  provided.  Appendix  2 is  a listing  of  the  main  pro- 


gram and  the  four  primary  subroutines.  Several  other  subroutines,  used 


only  for  plotting  purposes,  are  not  shown.  The  following  is  a list  of 


important  terms  along  with  some  discussion  of  each  term  as  necessary. 


These  are  presented  in  the  approximate  order  in  which  they  are 


encountered  in  the  program. 


Variable 


Function 


Main  Program 


IUN 

(array-card  17) 


This  defines  the  X values  of  the  velocity 
vs.  hori~ontal  distance  plots.  The 
location  of  these  X values  is  IUN*DEL 


PTIME 

(array-card  17) 


This  defines  the  times  at  which  plots 
will  be  drawn 


Read  statement  (card  19)  - first  data  card 


Initial  temperature 


POWER 


Power  of  laser  beam 


Beam  rcdius 


Size  of  time  steps 


Initial  pressure 


Number  of  components  in  the  gas  mixture 


Size  of  spatial  steps 


Length  of  side  of  cross-section  of 
the  square  enclosure 


Irradiation  time 


Read  statement  (card  38)  - next  N data  cards.  Each  subscript  M refers  to 
a different  component  in  the  gas  mixture.  Each  variable  is  an  array  of 
size  N. 


NAME 

RHOI 

ALPHA 

CONDI 

W 

CONC 

VNUE 

CP 

GRAV  (card  46) 

R (card  47) 

AVRHO  (card  49) 
AVALPH  (card  50) 

AVCON  (card  51) 

AVW  (card  52) 

A VNUE  (card  53) 
AVMU  (card  55) 

AVCP  (card  73) 

IMAX  (card  86) 

JMAX  (card  87) 

UN  (array-card  91) 
VN  (array-card  92) 


Name  of  component 

Initial  density  of  component  at 
TO  and  P0 

Absorption  coefficient  per  psi  and 
per  volume  percent  of  that  component 

Thermal  conductivity 

Molecular  weight 

Concentration  (mole  fraction) 

Kinematic  viscosity 

Specific  heat 

Acceleration  of  gravity  in  feet  per 
second^ 

3 

Universal  gas  constant  in  (psia-ft  / 
mole  °R) 

Gas  mixture  initial  density 

Gas  mixture  total  absorption  coefficient 
(c.r1)  at  the  initial  conditions 

Gas  mixture  thermal  conductivity 

Gas  mixture  molecular  weight 

Gas  mixture  kinematic  viscosity 

Gas  mixture  viscosity 

Gas  mixture  specific  heat 

Number  of  grid  spaces  or  elements 
in  vertical  (X)  direction 

Number  of  grid  spaces  or  elements 
in  horizontal  (Y)  direction 

Non-dimensional  local  vertical  velocity 
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Non-dimensional  local  horizontal  velocity 


RHON  (array-card  93) 
TN  (array-card  94) 

QN  (array-card  95) 

TF  (array-card  96) 

VORIP  (array-card  97) 

SI  (array-card  58) 
NPRINT  (card  100) 

PI  (card  102) 

PX  & PY  (card  107) 

KPTIME  (card  114) 

F (card  117) 

I TIME  (card  120) 
XCENT  (card  122) 

DTAUN  (card  131) 

DELN  (card  132) 

PRAND  (card  133) 
POWDEN  (card  144) 

GAM  (card  147) 


Non-dimensional  local  density 

Non-dimensional  local  temperature 
at  start  of  a time  step 

Non-dimensional  local  volumetric 
heating  rate 

Non-dimensional  local  temperature  at 
end  of  a time  step 

Non-dimensional  local  vorticity  at 
the  start  of  a time  step 

Non-dimensional  local  stream  function 

Whenever  the  number  of  time  steps  taken 
is  equal  to  NPRINT  (i.e.,  ITIME  .EQ. 
NPRINT)  then  the  pressure,  total  heating 
rate,  and  time  will  be  printed. 

Initial  pressure 

Used  to  establish  plotter  grid.  All 
points  in  the  plots  can  be  identified 
by  coordinates  (PX,  PY) . 

This  defines  the  subscript  value  for 
the  PTIME  array 

This  term  determines  the  size  of  the 
plots.  For  F=l,  the  plots  are  4"  x 8". 
The  size  of  plots  is  proportional  to  F. 

This  defines  the  number  of  time  steps 
preceding  the  current  one. 

This  determines  the  vertical  distance 
from  the  bottom  edge  of  the  enclosure 
to  the  center  of  the  enclosure. 

Non-dimensional  time  step  (At) 

Non-dimensional  spacial  step  (AX) 

Prandtl  number 

2 

Average  power  density  (watts /cm  ) based 
on  input  beam  power  and  radius. 

This  is  the  "dummy"  variable  used  in 
subroutine  MOTION  to  solve  the  equation 
of  motion  via  the  ADI  method. 


PN  (card  160) 


Pressure  at  the  end  of  each  time  step 


QTH  (card  175) 


TMAXPP  (card  206) 


The  actual  heating  rate  per  inch  of 
thickness  over  the  entire  beam  area. 
This  is  calculated  in  subroutine  HEAT 
for  one-half  of  the  beam  area. 

The  maximum  temperature  attained  at 
any  point  in  the  array.  The  program 
determines  this  value,  prints  it,  and 
automatically  plots  the  graphs  at  the 
next  time  step. 


I 

| 

I 

[ 


r 

R 

f 

I 

l 


Subroutine  HEAT 

This  subroutine  determines  which  grid  elements  are  to  be  heated 
(e.g.,  whether  they  be  in  beam  path)  and  calculates  the  local  heating 
rate  at  each  element.  The  proportional  intensity  of  the  beam  and  the 
coordinates  of  the  heated  edge  of  the  beam  are  also  determined  and  printed. 
The  total  volumetric  heating  rate  (QTH)  for  one-half  of  the  beam  and  per 
inch  thickness  is  also  determined.  Finally,  the  actual  power  and  actual 
beam  area  used  in  the  calculations,  which  differs  from  input  data  due  to 
beam  truncation,  and  numerical  approximation  errors,  is  determined  and 
printed.  Cards  20  and  21  are  changed  when  a "flat-top"  beam  is  considered. 


RAD  (card  15) 


ALPH  (card  18) 
AREA  (card  26) 


QPROP  (card  31) 


This  is  the  distance  from  an  element 
to  the  center  of  the  beam.  This  is 
used  to  determine  if,  and  how  much, 
heating  will  occur  at  that  element. 

Local  absorption  coefficient  per 
inch 

One-half  of  the  actual  cross-sectional 
area  in  the  beam  path.  Converted  to 
total  actual  beam  area  at  card  80. 

Ratio  of  the  intensity  of  the  beam 
at  any  point  to  its  intensity  at  the 
beam  center. 
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BX  & BY  (arrays-cards  62  & 6 ) These  define  the  heated  elements  which 

lie  on  the  edge  of  the  heated  (beam) 
area . 

QT  (card  81)  Actual  beam  power  used  in  program 

ACPOWD  (card  82)  Actual  power  density 


Subroutine  TEMP 


This  subroutine  solves  the  energy  equation  in  order  to  find  the  new 
temperature  array  TF. 


Subroutine  MOTION 

This  subroutine  solves  the  motion  equation  to  find  the  new  vorticity 
array  VORT.  This  new  VORT  array  is  used  to  find  the  SI  array  (stream 
function)  via  the  ADI  method  described  in  the  report.  Finally,  the  SI 
array  is  converted  to  local  velocity  arrays  UN  and  VN.  The  variable 
XCHEK  at  card  31  is  the  value  used  to  determine  if  suitable  convergence 
has  occurred  in  the  ADI  method  calculations. 

Subroutine  TRIDAG 

This  subroutine  solves  the  system  of  linear  equations  encountered 
in  the  ADI  method  ca] elation  of  stream  function  (SI)  from  vorticity 
(VORT)  in  subroutine  MOTION. 
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